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Abstract—The Ocellated Lizard, Timon lepidus (Daudin 1802) occupies the Mediterranean regions of 
southwestern Europe (Portugal, Spain, France, and the extreme northwest of Italy). Over the last decades, a 
marked decline in its population has been observed, particularly on the northern edge of its distribution. As a 
result, it is currently considered a threatened species, especially in France and Italy. In France, a national action 
plan for its conservation has been put in place. In this study, ecological niche modelling (ENM) was carried out 
over the entire area of France in order to evaluate the species’ potential distribution, more accurately define its 
ecological niche, guide future surveys, and inform land use planning so this species can be better taken into 
consideration. The modelling used data representing 2,757 observation points spread over the known range 
of the species, and 34 ecogeographical variables (climate, topography, and vegetation cover) were evaluated. 
After removing correlated variables, models were fitted with several combinations of variables using eight 
species distribution model (SDM) algorithms, and then their performance was assessed using three model 
accuracy metrics. Iterative trials changing the input variables were used to obtain the best model. The optimized 
model included nine determining variables. The results indicate the presence of this species is linked primarily 
to three climate variables: precipitation in the driest month, precipitation seasonality, and mean temperature 
in the driest quarter. The model was checked by a sample dataset that was not used to fit the model, and this 
validation dataset represented 25% of the overall field observations. Of the known occurrence locations kept 
aside to check the results, 94% fell within the presence area predicted by the modelled map with a presence 
probability greater than 0.7, and 90% fell within the area with a presence probability ranging from 0.8 to 1, which 
represents a very high predictive value. These results indicate that the models closely matched the observed 
distribution, suggesting a low impact of either geographical factors (barriers to dispersal), historical factors 
(dispersal process), or ecological factors (e.g., competition, trophic resources). The overlap between the 
predicted distribution and protected areas for this species reveals that less than 1% of the potential distribution 
area is protected by strong regulatory measures (e.g., national parks and natural reserves). The knowledge 
obtained in this study allows us to recommend some guidelines that would favor the conservation of this 
species. 


Résumeé.—Le lézard ocellé, Timon lepidus (Daudin 1802), occupe les régions méditerranéennes du sud-ouest 
de l’Europe (Portugal, Espagne, France, et extrême nord-ouest de l'Italie). Au cours des dernières décennies, 
un fort déclin des populations a été observé, particulièrement aux marges nord de sa distribution. Il est donc 
considéré comme une espèce menacée, spécialement en France et en Italie. En France, il bénéficie d’un 
plan national d’actions en faveur de sa préservation. La modélisation de sa distribution a été conduite sur 
ensemble du territoire national en vue d’estimer sa distribution potentielle, préciser sa niche écologique, 
orienter les prospections futures et permettre une meilleure prise en compte de l'espèce dans l'aménagement 
du territoire. Le travail de modélisation repose sur 2757 points d'observation répartis sur ensemble de la 
distribution connue de l'espèce, confrontés a 34 variables climatiques, topographiques, et de couvert végétal. 
Après suppression des variables autocorrélées, plusieurs combinaisons de variables ont été testées, et leur 
performances évaluées à partir de huit algorithmes SDM. Le meilleur modèle retient neuf variables, déterminées 
par l'algorithme ayant la meilleure performance. Les modèles montrent que la présence de espèce est 
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principalement déterminée par la sécheresse et la température estivale: précipitations au cours du mois le 
mois le plus sec, saisonnalité des précipitations et température moyenne des trois mois les plus chauds. La 
validation du modèle sur la base d’un échantillon totalisant 25 % du total des observations, non inclus dans le 
modèle, montre que 94 % des données de validation se placent dans laire potentielle au seuil de probabilité de 
0,7, et 90 % pour une probabilité comprise entre 0,8 et 1. Ceci donne une valeur prédictive tres élevée au modèle 
retenu. On constate une étroite concordance entre la distribution potentielle et la distribution réalisée, ce qui 
suggère une faible influence des facteurs géographiques (obstacles a la dispersion), historiques (processus 
de dispersion) ou écologiques (compétition, ressources trophiques, etc.). Le croisement cartographique entre 
Paire potentielle de espèce et les espaces protégés montre que moins de 1 % de laire potentielle est couverte 
par des mesures réglementaires fortes (parcs nationaux et réserves naturelles). En conclusion, le travail donne 
des orientations pour améliorer la connaissance de la distribution de espèce et des pistes de réflexion en 
faveur de sa conservation. 
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Introduction 


The success of conservation programs in protecting 
threatened species depends largely on the quality of 
the information related to the environmental condi- 
tions favorable to (or sought by) the species (Griffith et 
al. 1989; Souter et al. 2007; Fourcade et al. 2018). For 
this reason, ecological (or environmental) niche mod- 
els (ENMs) [Sillero 2011], also known as species dis- 
tribution models (SDMs) or habitat distribution models 
(HDMs), are increasingly used to inform conservation 
measures (Ferrier 2002; Graham et al. 2004; Araújo and 
Segurado 2004; Santos et al. 2009; Elith and Leathwick 
2009; Lyet et al. 2013; Jiang et al. 2014; Wan et al. 2016; 
Ianella et al. 2018). These models allow researchers to 
identify the factors which explain the distribution of a 
species (Austin et al. 1990; Vetaas 2002; Guisan and 
Hofer 2003; Jiang et al. 2014; Ferreira et al. 2013), to 
orient research toward zones where the species has not 
yet been identified (Engler 2004; Raxworthy et al. 2003; 
Lyet et al. 2013; GHRA-LPO Rhône-Alpes 2015; Ryberg 
et al. 2017), to identify the most favorable zones for the 
conservation of the species (Brito et al. 1996; Barbosa 
et al. 2003; Anderson and Martinez-Meyer 2004; Mufioz 
et al. 2005; Guisan et al. 2013; Lyet et al. 2013; Wang 
et al. 2016; Moradi et al. 2019; Sohrab et al. 2019), to 
evaluate the potential dispersion and gene flow between 
population patches (Guisan and Thuiller 2005), to model 
future changes in distribution (for example, based on cli- 
mate change) [Franklin 1998; Guisan and Hofer 2003; 
Araujo et al. 2006; Carvalho et al. 2011; Cheaib et al. 
2012; Ianella et al. 2018; Renwick et al. 2018], as well as 
to project the distribution into past scenarios (Sillero and 
Carretero 2013). 

Species distribution models also allow comparisons 
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between potential and actual distributions, two very use- 
ful concepts to consider in conservation biology. Put sim- 
ply (see Pulliam 2000 for a more detailed explanation), 
a potential niche corresponds to areas in which the cli- 
matic, terrain, and habitat conditions are theoretically fa- 
vorable to the species in the current conditions, whereas 
the realized niche takes into account historical and bi- 
otic factors that may explain the absence of the species 
within the ecological area defined by the fundamental 
niche. Comparing these two niches provides informa- 
tion about the historical processes that led to the cur- 
rent distribution, the dispersion capacity of the species, 
and the obstacles to its dispersion. This can illuminate 
the ecological factors that negatively influence the pres- 
ence of the species, such as the presence of competitors 
or predators, unsuitable habitats, insufficient trophic re- 
sources, a lack of host species, and others (see Guissan 
and Thuiller 2005). When used to model the distribution 
of a species in decline, SDMs can also provide informa- 
tion about the causes of decline, particularly by helping 
to differentiate the proportions due to global factors as 
opposed to regional or local factors (Jiang et al. 2014). In 
some cases, SDMs can even allow the extent of decline 
to be measured, by comparing the potential niche with 
the observed niche (Lyet et al. 2013; Ryberg et al. 2017). 

The Ocellated Lizard, Timon lepidus (Daudin 1802), 
is a good case study for this type of analysis. This species 
occupies the Mediterranean regions of southwestern Eu- 
rope (Portugal, Spain, France, and the extreme northwest 
of Italy) [Figs. 1 and 2]. At the edges of its distribution, 
it has faced a marked population decline over the last de- 
cades and is now considered a threatened species, espe- 
cially in France and Italy (Salvidio et al. 2004; Cheylan 
and Grillet 2005; Cheylan 2016). 

This species is closely linked to the Mediterranean 
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Distribution of Timon lepidus in France 


Fig. 1. Timon lepidus, adult male, Hérault, France. Photo by 
Jean Nicolas. 


climate and specific biotopes (Doré et al. 2015). More- 
over, due to its large size and thermal requirements for 
reproduction, it is particularly demanding regarding cli- 
matic conditions (Mateo 2011). The incubation period 
for its eggs is long, about 100 days, making use of the 
full period of warmer temperatures. The eggs hatch late, 
generally at the end of September or the beginning of 
October (Bischoff et al. 1984; Doré et al. 2015). Hence, 
two main factors drive the distribution of this species: its 
ecophysiology (thermal requirements linked to the indi- 
vidual’s bulk and to the incubation of the eggs) and its 
habitat requirements (dry habitats with little tree cover). 
In view of these needs, the Ocellated Lizard should theo- 
retically benefit from the warming temperatures recorded 
in Europe over the last 30 years (see Prodon et al. 2017). 
Yet observations show a rather widespread decline of this 
species, particularly on the northern edge of its distribu- 
tion range, which is inconsistent with the expected situa- 
tion in warming conditions (Salvidio et al. 2004; Cheylan 
and Grillet 2005; Doré et al. 2015). This raises questions 
regarding the causes of the decline of this species and, at 
first glance, suggests a hypothesis that local effects pre- 
dominate over global effects. 

The use of SDMs allows the study of interesting bio- 
geographical questions about this species. Native to the 
Iberian Peninsula, the Ocellated Lizard colonized France 
and the extreme west of Italy along the Mediterranean 
coast (Sillero et al. 2014; Doré et al. 2015). This colo- 
nization involved overcoming considerable physical ob- 
stacles, including rivers, mountains, and forests. Niche 
modelling can provide information about constraints that 
limit dispersion; that is, whether the current distribution 
boundaries of the species are of a climatic nature (thus 
ecophysiological) or a physical nature (due to obstacles 
to dispersion). The same question applies to the species’ 
colonization of the French Atlantic coast: Are the isolat- 
ed populations along this coast the result of a process of 
decline linked to the progressive degradation of habitat 
or to climatic constraints? The responses to these ques- 
tions can be found by comparing the expected distribu- 
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Fig. 2. Timon lepidus, juvenile, Hérault, France. Photo by Jean 
Nicolas. 


tion with the observed distribution; that is, by comparing 
the potential niche with the realized niche. 

This study used species distribution modelling to in- 
vestigate the following questions: (1) Does the observed 
distribution of the Ocellated Lizard match its potential 
distribution? If not, why not? (2) Which variables best 
explain the distribution of this species: climate, terrain, 
land use, or other factors? (3) Why is this species retreat- 
ing at the edges of its distribution range, in contrast to 
what might be expected based on climatic changes? (4) 
Are the distribution boundaries of this species conditional 
on either climate, physical barriers, or ecological causes? 
(5) Which zones are potentially the most favorable for 
the conservation of this species? (6) Where should future 
surveys be carried out to improve our understanding of 
the distribution range? (7) Based on these findings, what 
conservation strategy should be implemented for the 
conservation of the species? 


Materials and Methods 
Data 


Ocellated Lizard dataset. This study used observations 
(presence data only, Brotons et al. 2004) from an exhaus- 
tive database that includes most of the occurrences ob- 
served in France between 1970 and 2016. As the objec- 
tive was essentially practical, i.e., to identify the areas 
of the potential presence of the species with the aim of 
its conservation, this study did not consider taking into 
account the entire distribution of the species as relevant 
to building the model. This would have led to further 
complications, such as the need to consider markedly di- 
vergent genetic lines and, as a result, ecophysiological 
adaptations specific to the regions that host these genetic 
lines. 

The data were collected by a number of organizations 
and individuals over a period of 46 years. Before inte- 
grating the data into the models, the records were verified 
and cross-validated, keeping only precisely georefer- 
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Fig. 3. Localization of the total presence data for the Ocellated Lizard TOn merre collected i in Benes during the period 1970 to 


2016. 


enced locations (either data captured by GPS or observa- 
tions recorded with a spatial positioning error of less than 
100 m). This resulted in presence data for a total of 5,521 
locations spread over southern France. From this, a sam- 
ple dataset was extracted for modelling purposes. First, 
only occurrences for the period that corresponded to the 
vegetation variable used in the model were selected, tak- 
ing into account the development and continuity of the 
vegetation cover over the years. To fit this requirement, 
the presence data were narrowed down to a sourcing pe- 
riod of 16 years, which included 4,282 observations from 
2000 to 2016. The distribution of occurrence data for this 
period covers the whole area of study, though data prior 
to the year 2000 were not included (Fig. 3). 

Secondly, the dataset was filtered to avoid spatial bias 
due to oversampling at particular locations, as field in- 
vestigations conducted for environmental impact assess- 
ments and other monitoring programs led to a higher con- 
centration of data at certain sites. Therefore, the density 
of points per km? was evaluated to identify zones subject 
to sampling bias (Fig. 1 in Supplementary Materials), 
using kernel density calculation. Through this analysis, 
zones where the point density ranged from 2—50 points 
per km? were determined, and for these zones, a single 
record was retained per km”. After filtering the data by 
density, a total of 2,757 occurrences were retained, pro- 
viding presence data that was well distributed over the 
study area (Fig. 4). Finally, of these 2,757 valid records, 
75% were randomly selected for modelling, with the re- 
maining 25% serving as an independent source to check 
(validate) the results (Hirzel and Guisan 2002) [Fig. 4]. 
During the modelling procedure, the database consisting 
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of 75% of the occurrences was itself divided into training 
data (80%) and testing data (20%). 


Environmental data. A set of ecogeographical variables 
was used to model species distribution, taking into ac- 
count the ecological requirements of the species (Guisan 
and Thuiller 2005). A group of 34 variables was evalu- 
ated through several iterations in order to identify and 
include the most relevant variables (Table 1). 

As a first step, environmental variables were obtained 
from available sources at the appropriate spatial and the- 
matic resolutions. This study used the CHELSA data- 
base, version 1.1 (Climatologies at High Resolution for 
the Earth’s Land Surface Areas, Karger et al. 2017). This 
database includes a set of bioclimatic variables, with 
monthly mean temperature and precipitation patterns, 
for the time period 1979—2013, which corresponds to the 
time range of the species occurrence data. The CHEL- 
SA database is an alternative source to the widely used 
WorldClim global climate database, as both derive their 
bioclimatic variables from the monthly minimum, maxi- 
mum, and mean temperatures, as well as precipitation 
values. However, as described by Karger et al. (2017), 
the CHELSA variables include additional corrections, 
such as monthly mean and station bias, wind effect and 
valley exposition, as well as correction for orographic 
effects on precipitation. Additionally, as CHELSA is a 
recently released product, this study allowed evalua- 
tion of its potential for species distribution modelling. 
In a recent study, Karger et al. (2017) highlighted differ- 
ences observed at a large scale between WorldClim and 
CHELSA models, with the latter leading to a significant 


December 2019 | Volume 13 | Number 2 | e213 


Distribution of Timon lepidus in France 


3 few 
e Ta on 
e s 
. = 
a 
l rou B 
H Gace 
} S taa 
l [z eA 
aad ° 3 Se 
pm Man A A 
c Pale 22:20 


“ome 


Andotra 


0 2550 100km 
i > 


Sources Esn, 
METI Esri China (Hang Kong), swisstopo, Mapmyincia, 


75% of occurence data used 
in the model 


25% of occurence data used 
for validation 


Esn, HERE Detorme, ntermap, increment P Ci 
, Swiss o, Mapmyincia, q 


Fig. 4. Localization of occurrence data for the Ocellated Lizard Timon lepidus, showing the points in the dataset used for the model 


(red) and those used to check the model (blue). 

improvement in SDM performance. The CHELSA bio- 
clim data grid has a 30-arc-second pixel resolution, cor- 
responding to a spatial resolution of 0.00833 decimal de- 
grees at the equator. Applied to the study area in southern 
France, the bioclimatic variables have a spatial resolution 
of less than 1 km’, with each grid roughly covering an 
area of 700 x 900 m. 

As the CHELSA database does not include variables 
reflecting solar radiation and its impact on climatic hu- 
midity or aridity, data produced by the Consultative 
Group on International Agricultural Research (CGIAR) 
were used (Fick and Hijmans 2017). Data available at the 
same resolution as the CHELSA dataset included mean 
annual solar radiation, potential evapotranspiration, and 
a global aridity index calculated from the ratio of mean 
annual precipitation to mean annual potential evapo- 
transpiration. 

During the process of selecting and consolidating the 
ecogeographical variables, another climate data source 
was used to compare the results and identify the best input 
that would optimize performance. This climate dataset 
(Joly et al. 2010) is based on data from Météo France’s 
weather stations and gathers a number of indicators that 
are not provided in the CHELSA dataset, such as the 
number of days with a temperature above 30°C or below 
-5 °C, as well as other parameters related to temperature 
ranges and seasonal variations over the year. 

To increase model performance, topographic data 
were also included. The elevation was obtained from the 
EU-DEM v1.1 dataset, produced in the framework of 
the European Commission’s Copernicus program (Bash- 
field and Keim 2011), which includes a digital elevation 
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model captured in 2011 and projected in ETRS89-LAEA 
(EPSG code 3035), with a spatial resolution of 25 m. 
This allowed the modelling to take altitude into account, 
as well as slope and aspect, calculated in projected coor- 
dinate systems (all in meters). 

A variable representing vegetation cover is also re- 
quired, in order to identify natural habitats suitable for 
the Ocellated Lizard. The vegetation cover indicator se- 
lected was the normalized difference vegetation index 
(NDVI) generated from the MODIS (Moderate resolu- 
tion Imaging Spectro-radiometer) Terra satellite sensors 
(Huete et al. 2010), which has a spatial resolution of 250 
m. Data captured in the first week of July were used, as 
this corresponds to the optimal season for identifying 
the permanent vegetation cover of interest in this study. 
Vegetation indices in summer can highlight chlorophyll 
produced by permanent vegetation, whereas vegetation 
indices calculated in spring are influenced by the annual 
growth of the herbaceous stratum. 

Considering the possibility that land cover has 
changed over time, the importance of changes in veg- 
etation cover over the years was evaluated. The NDVI 
values were compared for each year from 2000 to 2016, 
calculating the variance and the standard deviation of 
the mean NDVI value for this 16-year period (Fig. 5). 
The results indicated that the change in permanent veg- 
etation cover in the area of interest has been negligible, 
especially in areas where the Ocellated Lizard has been 
observed. Over 16 years, for the month of July, the stan- 
dard deviation of the yearly NDVI values relative to the 
NDVI mean value was below 0.075 for 98.9% of the en- 
tire studied area, below 0.05 for 88.9% of the area, and 
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Table 1. Description of the 34 variables evaluated. 


Variable Description 
Biol Annual mean temperature 
Bio3 Isothermality (BIO2/BIO7) (* 100) 
Bio4 Temperature seasonality (standard deviation *100) 
Bio5 Max temperature in warmest month 
Bio8 Mean temperature in wettest quarter 
Bio9 Mean temperature in driest quarter 
Biol0 Mean temperature in warmest quarter 
Biol3 Precipitation in wettest month 
Biol4 Precipitation in driest month 
Biol5 Precipitation seasonality (coefficient of variation) 
Biol6 Precipitation in wettest quarter 
Biol7 Precipitation in driest quarter 
Bio 18 Precipitation in warmest quarter 
Biol9 Precipitation in coldest quarter 
Alt Elevation 
Slope Slope 
Aspect Slope orientation 
NDVI Normalized vegetation index 
TCD Tree cover density 
srad Mean annual solar radiation 
PET Global potential evapotranspiration 
Global aridity index (mean annual precipitation / mean 
Aridity annual PET) 
TMO Annual mean temperature 
TMN Number of days with temperature below -5 °C 
TMX Number of days with temperature above 30 °C 
TAM Annual temperature range 
DEH Inter-annual temperature variability in January 
PEE Inter-annual temperature variability in July 


Source Period Resolution 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~] km? 
CHELSA v1.1 1979-2013 ~l km? 
CHELSA v1.1 1979-2013 ~l] km? 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~l km? 
CHELSA v1.1 1979-2013 ~] km? 
CHELSA v1.1 1979-2013 ~1 km? 
CHELSA v1.1 1979-2013 ~1 km? 
EU-DEM v1.0 European 
Commission Copernicus program 2011 25m 
Calculated from EU-DEM v1.0 2011 25m 
Calculated from EU-DEM v1.0 2011 25m 
MODIS TERRA (NASA) July 2012 250 m 
Copernicus 2012 20m 
WORLDCLIM v.2 1970-2000 ~1 km? 
CGIAR 1950-2000 ~] km? 
CGIAR 1950-2000 ~1 km? 


ThéMA - CNRS-UMR 6049 1971—2000 250 m 


ThéMA - CNRS-UMR 6049 1971—2000 250 m 


ThéMA - CNRS-UMR 6049 1971—2000 250 m 
ThéMA - CNRS-UMR 6049 1971—2000 250 m 


ThéMA - CNRS-UMR 6049 1971—2000 250 m 


ThéMA - CNRS-UMR 6049 1971—2000 250 m 


Variance between January precipitation and monthly mean 


PDH precipitation 

Variance between July precipitation and monthly mean 
PDE precipitation 
PJH Number of rainy days in January 
PEH Inter-annual precipitation variability in January 
PEE Inter-annual precipitation variability in July 


ThéMA - CNRS-UMR 6049 1971-2000 250 m 


ThéMA - CNRS-UMR 6049 1971-2000 250 m 
ThéMA - CNRS-UMR 6049 1971-2000 250 m 


ThéMA - CNRS-UMR 6049 1971—2000 250 m 


ThéMA - CNRS-UMR 6049 1971—2000 250 m 


Variation between autumn (September and October) and 


PRA July precipitation 


only reached a maximum of 0.2 for 1.1% of the area. An 
assessment of the changes in vegetation cover within the 
study area allowed the selection of the NDVI mean value 
for the 16-year NDVI dataset (from 2000 to 2016). 

All the selected variables were aggregated to a spa- 
tial resolution of ~1 km, using bilinear resampling tech- 
niques. 
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ThéMA - CNRS-UMR 6049 1971-2000 250 m 
Methodology 


The performance of models was explored based on dif- 
ferent combinations of ecogeographical variables using 
an iterative approach (Heikkinen et al. 2006). To do this, 
a model was first fitted with a selected set of variables 
(Bucklin et al. 2014) and its performance was measured. 
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Table 2. Input variables for the six models. See Table 1 for descriptions of variables. 


Code Model description 
M1 six selected bioclimatic variables 
M2 eight random bioclimatic variables (blind test) 
M3 12 climatic variables from Météo France 
M4 M1 variables + solar radiation + PET + aridity 
M5 M1 variables + altitude, slope, aspect 
M6 M1 variables + altitude, aspect + vegetation 


The inputs were then optimized by testing the model with 
other sets of variables, using a stepwise procedure to in- 
clude or exclude variables one by one. Every iteration 
was evaluated and outcomes associated with the chosen 
parameters were recorded. This resulted in six models 
based on different combinations of variables (Table 2). 
While a comparison of bioclimatic data sources in 
terms of usefulness for species distribution modelling is 
beyond the scope of this study, a deliberate choice was 
made to use various sources in the model; specifically, 
bioclimatic data from entirely different inputs. This itera- 
tive and multiple-source approach increases the chances 
of obtaining a successful model by allowing the selection 
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for the month of July was calculated from MODIS data captured on the first 


week of July every year from 2000 to 2016. 
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Variables 
Bio4, Bio9, Biol0, Biol4, Biol5, Biol6 
Bio3, Bio4, Bio5, Bio8, Bio9, Biol3, Biol7, Biol9 


TMO, TMN, TMX, TAM, TEH, TEE, PDH, PDE, PJH, 
PEH, PEE, PRA 


[Biol, Bio4, Bio9Bio15, Biol6] + srad + PET + Aridity 


[Bio4, Bio9, Biol0, Biol4, Biol5, Biol6] + Alt, Slope, 
Aspect 
[Bio4, Bio9, Biol0, Biol4, Biol5, Biol6] + [Alt, Aspect] 
+ NDVI 


of the most accurate and relevant variables. It also helps 
to validate overall model performance, by providing keys 
for analyzing the suitability of the model in terms of bio- 
climatic variables favorable to the species. 


Model descriptions 


Model I: Variables selected from the CHELSA climate 
dataset based on ecological assumptions. The first 
modelling trial used a set of ecogeographical variables 
identified based on expert knowledge of the biology 
and behavior of the species. As the Ocellated Lizard is a 
Mediterranean species that prefers long periods of warm 
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Fig. 5. Change in vegetation cover from 2000 to 2016: standard deviation of the yearly NDVI values relative to the NDVI mean 


value for the month of July. 
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Table 3. Variable correlation values for Model 1. 


Jorcin et al. 


bio4 R bio9_R biol0_R biol4 R biol5_R biol6 R 
bio4 R 1.00000000 -0.07440553 0.3081616 -0.2788248 0.1672821 -0.2605788 
bio9_R -0.07440553 1.00000000 0.5115825 -0.4025089 0.3664137 -0.0604624 
biolO_R 0.30816159 0.51158248 1.0000000 -0.6782349 0.6210210 -0.1046921 
biol4_ R -0.27882481 -0.40250891 -0.6782349 1.0000000 -0.6896546 0.2434781 
biol5_R 0.16728213 0.36641372 0.6210210 -0.6896546 1.0000000 0.2159100 
biol6 R -0.26057884 -0.06046240 -0.1046921 0.2434781 0.2159100 1.0000000 


temperatures (for ecological reasons linked to its repro- 
ductive cycle), this model considered the bioclimatic 
variables that would best represent this need. Variables 
with a recognized influence on a species’ ecology are 
generally expected to lead to more accurate predictions 
in SDMs (Austin 2002). This initial SDM incorporated 
the following variables from the CHELSA database: 
Bio4, Bio9, Biol0, Biol4, Biol5, and Biol6 (Table 3). 
The Biol6 variable was selected as it represents precipi- 
tation in the wettest quarter, so it would be a good in- 
dicator of the aridity of the environment. This variable 
differentiates regions according to precipitation patterns 
by indicating rainfall occurring during the wettest period 
of the year, thus allowing locations with lower rainfall 
to be identified. After this a priori selection, the correla- 
tions between variables were measured (Table 3). As the 
variables were not correlated, this choice was validated. 


Model 2: Non-correlated variables selected from the 
overall CHELSA climate dataset. For the next model, the 
correlations between the 19 variables obtained from the 
original CHELSA database for all presence-data loca- 
tions were calculated, removing one of each pair of high- 
ly correlated variables (those with a correlation coeffi- 
cient greater than 0.75) [Doorman 2012]. As 11 of the 19 
input variables were correlated, only the remaining eight 
were retained, with no consideration of any presumed 
ecological significance. Thus, this trial was considered 
a blind test, run on a statistical basis only rather than on 
prior knowledge of the input variables. Problems of col- 
linearity between variables were identified and dealt with 
using variance inflation factors (VIF) with the R usdm 
package (Naimi et al. 2014). A VIF was calculated for 
each explanatory variable, and those with a VIF greater 
than 10 were removed. The correlation coefficients of 
the remaining variables ranged between -0.015 and 0.75 
(Table 1 in Supplementary Materials). 


Model 3: Non-correlated variables from the Météo 
France climate dataset. As an alternative to the CHEL- 
SA climate database, this model used a climate dataset 
obtained from Météo France weather stations (Joly et al. 
2010). Correlation tests showed that out of 14 variables 
from the Météo France dataset, only two were highly 
correlated, with a correlation coefficient greater than 
0.75. Therefore, the 12 non-correlated variables were in- 
cluded in the model, with each having a potential effect 
on model performance. 
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Model 4: The variables for Model 1 with the addition 
of solar radiation, evapotranspiration, and aridity. 
This model included the six climate variables selected 
for Model 1 in addition to climatic parameters that are 
potentially important to the ecology of the species. To 
reflect discriminating factors related to the Mediterra- 
nean climate, this model used mean annual solar radia- 
tion (obtained from WorldClim), as well as the Global 
Potential Evapotranspiration and Global Aridity Index 
(obtained from CGIAR) [Fick and Hijmans 2017]. The 
Global Aridity Index consists of mean annual precipita- 
tion divided by mean annual potential evapotranspiration 
(Zomer et al. 2008). An assessment of whether the per- 
formance of the model increased with these additional 
variables was made. 


Model 5: The variables for Model I with the addition 
of topographic variables. This model included the six 
climate variables selected for Model 1 along with three 
additional topographic variables: elevation, slope, and 
aspect (orientation). Topographic parameters were ex- 
pected to improve model performance (Humboldt and 
Bonpland 1805). 


Model 6: The variables for Model I with the addition 
of selected topographic variables and a vegetation 
variable. The final model included the six climate vari- 
ables selected for Model 1 along with two additional 
topographic variables (elevation and aspect) as well as 
NDVI. The two topographic parameters were retained 
because of the gain in performance obtained by add- 
ing them to the bioclimatic parameters. The addition of 
the NDVI helped to account for vegetation cover as a 
contributing variable in the model, as this is a valuable 
parameter in identifying the natural habitat of the spe- 
cies. Assuming the distribution of the Ocellated Lizard 
is linked to this species’ preferences in terms of land 
cover, vegetation density is expected to help differenti- 
ate areas of occurrence from areas of absence (Wilson 
et al. 2013). 


SDM methods 


To maximize SDM accuracy, all six models were run 
with eight statistical algorithms (Bucklin et al. 2014), 
regression-based machine learning, and classification 
methods. With each algorithm resulting in different pre- 
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Table 4. Algorithms used in the species distribution modelling. 


Code Description 

BRT Boosted Regression Trees 

CART aes ms Regression Trees for 
GAM Generalized Additive Model 

GLM Generalized Linear Model 

MARS Multivariate Adaptive Regression Spline 
MAXLIKE Maximum Likelihood 

RF Random Forests 

SVM Support Vector Machine 


dictions, the objective was to identify the method that 
achieved the best accuracy (Elith et al. 2006). Testing 
eight algorithms also allowed the evaluation of the over- 
all modelling approach (Table 4). As well as analyzing 
discrepancies between models in terms of performance, 
model congruence was examined to consolidate the 
conceptual approach (Li and Wang 2012). The model- 
ling methods used were Generalized Linear Modelling 
(GLM, Guisan and Zimmerman 2000), Generalized Ad- 
ditive Modelling (GAM, Guisan and Zimmerman 2000), 
Multivariate Adaptive Regression Spline (MARS, Elith 
and Leathwick 2007), Maximum Likelihood (MAX- 
LIKE), Classification and Regression Trees for Machine 
Learning (CART, Breiman et al. 1984), Boosted Regres- 
sion Trees (BRT, Elith et al. 2008), Support Vector Ma- 
chine (SVM, Drake et al. 2006), and Random Forests 
(RF, Breiman 2001). 

The presence-absence models were used with the ob- 
jective of predicting the presence probability of the Ocel- 
lated Lizard and mapping its distribution accordingly. 
Lacking absence data, pseudo-absences were used to 
run the models. Pseudo-absence data was generated with 
the R sdm package (Naimi and Araújo 2016), which has 
the advantage of providing a pseudo-absence selection 
process calibrated to SDM performance by considering 
presence data. Other studies have found that randomly 
selected pseudo-absences yield the most reliable models 
(Barbet-Massin et al. 2012). The models were fit by as- 
signing a number of pseudo-absences weighted to pres- 
ences (Barbet-Massin et al. 2012), with an equal number 
of presences and absences. 

Each type of SDM methodology has particular strengths 
and limitations in the way it accommodates the responses 
to predictors, as well as how it deals with missing observa- 
tions. For example, linear regression fits linear functions 
relating a response variable to one or more predictor vari- 
ables, where this relationship can be approximated by a 
straight line (Ferrier et al. 2002), whereas machine learn- 
ing offers more complex classification algorithms that 
accommodate non-linear variable interactions (Salas et 
al. 2017). All algorithms were tested with replicated sub- 
sampling of 20% of the occurrence dataset. 
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Model evaluation methods 


The contribution of SDMs in understanding the 
geographical distribution and abundance of a species 
depends on the level of reliability offered by the model 
(Barry and Elith 2006). The prediction accuracy must 
be assessed to determine the model’s suitability (Liu 
et al. 2009). Models can be judged on their capacity to 
discriminate presence from absence, which is measured 
by the number of false positive and false negative 
predictions. Several statistical indicators can be used as 
metrics to evaluate model performance (Fielding and 
Bell 1997). To assess the results here, the area under the 
Receiver Operating Characteristic (ROC) curve (AUC) 
value was used, as well as the correlation coefficient 
(COR) and the True Skill Statistics (TSS) value (Bradley 
1997). The AUC value provides a single measure of 
model performance, showing the model’s ability to 
rank a randomly chosen presence observation higher 
than a randomly chosen absence observation (Liu et 
al. 2009). These values can range between 0 and 1; and 
models producing AUC values of 0.75 are regarded as 
reliable, 0.8 as good, and 0.9 to 1 as having excellent 
discriminating ability (Franklin 2009). The TSS is 
presented as an improved measure of model accuracy, 
defining the average of the net prediction success rates 
for presence sites and for absence sites (Allouche et al. 
2006). The COR value allows another performance index 
comparison between models, and helps to validate the 
results obtained by the AUC and TSS methods (Elith et 
al. 2006). This study also ran null models, which make 
predictions in the absence of a particular ecological 
mechanism (Harvey et al. 2003), to assess the random 
probability hypothesis. The AUC values from the null 
models ranged from 0.51 to 0.53, thus corresponding to 
what would be expected by chance (Raes and ter Steege 
2007). The distribution generated with a null model 
significantly differed from the other modelling results, 
with a predicted presence spread smoothly over most of 
the study area (Fig. 2 in Supplementary Materials). 

As well as evaluating model performance with 
statistical indicators, the results were checked against 
the validation dataset (Anderson et al. 2003), consisting 
of the 25% of available presence data set aside for 
validation. The data distribution showed that 76% 
of the observation records within the presence range 
correspond to a threshold value of 0.9 to 1, and 86% were 
above a threshold of 0.8 (Table 5). Of the 690 occurrence 
locations set aside to check the results, 92% fell within 
the area predicted by the modelled map, in line with a 
presence probability threshold of 0.70. 

Predictive maps generated by SDMs provide the oc- 
currence probability of the species on a 0 to 1 scale. 
Threshold determination is a key step in transforming in- 
dices of suitability to binary predictions of species pres- 
ence or absence (Nenzen and Araujo 2011). Threshold 
definition can be subjective or objective (Manel et al. 
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Table 5. Correspondence between observations and presence 
probability for the validation dataset. Values above a threshold 
of 0.7 are shown in bold. 


Probability threshold eae He 
0.9 to 1 2,219 80.4 

0.8 to 0.9 266 9.7 

0.7 to 0.8 114 4.1 

0.6 to 0.7 80 a 

0.5 to 0.6 43 ie 

0.4 to 0.5 13 oe 

0 to 0.4 22 ae 


1999), and in many methods, the point at which sensitiv- 
ity (true positive rate) and specificity (true negative rate) 
are equal can be chosen to determine the threshold. In 
this case, a value of 0.7 was chosen. 


Results 


All SDMs for the differing sets of variables performed 
well, demonstrating high predictive power (Table 6). The 
overall mean AUC value was 0.86, the mean TSS was 
0.67, and the mean COR was 0.61. These results con- 
firmed the hypothesis that SDMs based on bioclimatic 
variables could provide valuable results for the Ocel- 
lated Lizard. Moreover, the agreement between results 
validates the overall modelling methodology, including 
the quality of the sampling dataset and the geographical 
extent of the study. Overall, the initial trial made with 
Model 1 showed that it was a satisfactory model, with 
an AUC value of 0.91, a TSS value of 0.72, and a COR 
value of 0.78. The subsequent iterative trials carried out 


with different combinations of variables further increased 
model accuracy, validating the most useful variables and 
helping to rank their contributions. 

A comparison of the overall results led to the selection 
of Model 6 generated with the RF method as the most 
accurate model. The predictors of this model, based on 
selected bioclimatic variables with additional topograph- 
ic and vegetation parameters, achieved the best perfor- 
mance using virtually all modelling methods, with maxi- 
mum performance obtained from the RF method. This 
combination of variables and modelling method resulted 
inan AUC value of 0.98, a TSS value of 0.85, anda COR 
value of 0.88. An analysis of variance (ANOVA) was 
performed to confirm the validity of the chosen combina- 
tion of variables. ANOVA results on Model 6 showed the 
lowest p-values for all variables, with a minimum value 
of 2.74e-14 and a maximum value of 0.020393 (Table 7). 
In contrast, models 1 to 5 each included a variable with 
a p-value > 0.05. As a study by Wood et al. (2016) men- 
tioned that the Akaike information criterion (AIC) works 
reasonably well for model selection, the AIC between 
models were also compared here, and this comparison 
showed that Model 6 had the lowest AIC values (Table 2 
in Supplementary Materials). 


Contribution of Variables 


The contributions of each variable in Model 6 were 
ranked, identifying four critical parameters, as well as 
one secondary factor (Table 8). The four main factors 
that most influenced model performance were (according 
to their relative importance): precipitation in the driest 
month, temperature seasonality, mean temperature in the 
driest quarter, and NDVI. Precipitation seasonality was 


Table 6. Comparison of model results based on different modelling methods and assessments of model performance. The values in 
bold, for Model 6 and RF, indicate the results with the highest accuracy. 
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Fig. 6. Predictive modelling maps showing presence probability of the Ocellated Lizard (Timon lepidus) in the study area. 


also a contributing variable, but only to a minor extent. 
Climatic factors associated with dry and warm weather 
conditions seem to play a determining role in the spatial 
distribution of the Ocellated Lizard, with precipitation in 
the driest month being the primary contributing factor. 
The contributions of these variables can be interpreted as 
a reflection of the species’ ecological needs, especially as 
related to its reproductive cycle. 


Predictive Habitat Suitability Maps 


The predictive maps (Fig. 6) show very slight differences 
between the six models. Oléron Island, where the most 
northerly currently known population of this species is 
found, is included in all models, but with different prob- 
ability ranges. Its presence is predicted on the whole is- 
land in models 1, 2, 4, and 5, but only on part of the island 
in models 3 and 6. All models show a clear link between 
the ‘Mediterranean population’ and the ‘Lot population’ 
(in a region lying northwest of the Mediterranean), with 
minor variations in the continuity of the species distri- 
bution between these two ‘populations.’ In the different 
maps, the penetration of the species into the Rhône Val- 
ley appears more or less extended, and the fragmentation 
of the “Lot population’ is more or less pronounced. Apart 
from these details, the maps based on the six models are 
extremely consistent. 
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On the map generated by Model 6 (Fig. 7), most of 
the locations historically occupied by the species (i.e., 
for which evidence exists of its disappearance) are along 
the Atlantic coast (in the Nouvelle-Aquitaine region), the 
area where the distribution of the species is the most lim- 
ited and fragmented (Fig. 8). The disappearances of the 
three Mediterranean populations correspond to very spe- 
cific cases: two islands (Ratonneau and Porquerolles is- 
lands), where the species is likely to have disappeared as 
a result of the introduction of predators (Cheylan 2016), 
and a population in the Rhône delta in the Camargue, 
where the decline of rabbits has transformed the envi- 
ronment, leading to the disappearance of the Ocellated 
Lizard in this area (Doré et al. 2015). 


Identification of Knowledge Gaps 


The potential distribution predicted by Model 6 indi- 
cates that knowledge gaps regarding the observed pres- 
ence of this species are considerable, not only in the 
Mediterranean distribution range, but also in its pe- 
riphery (Fig. 9). Allowing a buffer zone with a 5-km 
radius around each observation location, only 74% of 
the area of predicted presence (based on Model 6 with 
a presence probability threshold of 0.70) is confirmed 
by actual observations; observation data is lacking for 
26% of the area. In the core of the distribution range, 
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Fig. 7. The predicted distribution map generated by Model (in pink) with specific locations of presence data (from observations) 


from the dataset (black dots). 


information is missing in several areas of the regions 
of Provence and in a few areas of Languedoc-Roussil- 
lon. This is particularly the case for several noteworthy 
zones. 


e Var department (e.g., around Fayence, Draguig- 
nan, Bargemon, Seillons-source-d’Argens, and 
Saint-Maximin-la-Sainte-Baume) 

e Southeast of Alpes-de-Haute-Provence (e.g., the 
Valensole plateau, the Asse-Puimichel valley, 
and the lower Bléone valley) 

e Southeast and northwest of the Vaucluse depart- 
ment (e.g., Pertuis, La Tour-d’Aigues, Carpen- 
tras, Uchaux, Sainte-Cécile-les-Vignes, and 
Valréas) 

e Drôme department (Valence plain, Monteéli- 
mar plain, Tricastin, and Baronnies Proven- 
cales) 

e Gard department (e.g., Saint-Quentin-la-Poterie, 
Bagnols-sur-Ceéze, Alès, and La-Grand-Combe) 

e Southeast of Ardèche (e.g., Saint-Thomé, Vil- 
leuneuve-de-Berg, and Saint-Marcel-d’ Ardèche) 

e South and far northwest of the Lozère (e.g., 
Ayssèries, Faveyrolles, La-Bastide-Solages, Sé- 
brassac, Decazeville, and Claunhac) 


In western Languedoc-Roussillon, it would seem rel- 
evant to look for the Ocellated Lizard in several zones 
of the Aude department (e.g., Carcassonne, Labécède- 
Lauragais, Rouffiac-des-Corbières, and Palairac), 
and in several areas of the Pyrénées-Orientales (e.g., 
Saint-Paul-de-Fenouillet, and La Trinité) and the 
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Hérault (e.g., Pézenas and Aigues-Vives). Outside of 
the Mediterranean region, areas where this species 
would merit further survey efforts are more numerous, 
notably in the Aveyron, Tarn, Haute-Garonne, Tarn-et- 
Garonne, Lot, and in Dordogne, where the potential 
distribution area is very fragmented and the natural 
habitats small (Berroneau 2012; Pottier et al. 2017). 
To address this, surveys could be carried out in several 
areas of the Tarn and the Lot (e.g., Mazamet, Roque- 
courbe, Larroque, Lacapelle-Marival, Prayssac, Gour- 
don, and Martel), in southwest Corrèze (e.g., Tulle, 
Taurisson, and Saint-Aulaire), in southeast Dordogne 
(e.g., Carsac-Aillac, Hautefort, and Rouffignac-Saint- 
Cernin-de-Reilhac), and in eastern Tarn-et-Garonne 
(e.g., Caylus). All of these areas possess a high pres- 
ence probability of this species according to the model 
results, so it is likely that inadequate surveying ex- 


Table 7. ANOVA produced by the Generalized Additive Model 
for Model 6. Approximate significance of smooth terms. 


Edf Ref.df Chi.sq p-value 
s(bio4_R) 7.928 8.263 84.87  2.74e-14 
s(bio9_R) 7.445 T853 38.47  5.88e-06 

s(biol0_R) 6.469 6.779 42.75  7.371e-07 
s(biol4_R) 7.608 8.248 18.61 0.020393 
s(biol5_R) 6.071 7:297 18.72 0.008308 
s(biol6_R) 7.944 8.688 18.79 0.018410 
s(alt250_R) 6.279 6.819 34.27 9.87e-06 
s(orient250 R) 3.392 4.24] 20.47 0.000507 
s(median MO- 2.734 3.462 40.17 3.36e-08 
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Fig. 8. Location of populations that have disappeared in relation to the predicted presence map generated by Model 6 (blue dots). 


plains the current absence of data in these locations. 


Contribution of Protected Areas to the Conservation 
of the Species 


A comparison of the potential niche (Model 6 with a 
presence probability of > 0.70) and the main national 
protected areas allowed an assessment of the contribu- 
tion of the current network of natural reserves to the con- 
servation of the species (Table 9). 

The map in Fig. 10 shows the protected areas that 
contribute most to the conservation of the species are the 
regional nature parks (representing 22% of its predicted 
niche) and the two Natura 2000 zones (15% and 16%). 
The other types of protected areas contribute to a much 
more limited extent, covering between only 0.1% and 2% 
of the potential distribution area of 35,805 km’. As these 
protected areas often overlap (regional nature parks of- 
ten include nature reserves, and are generally also part of 
the Natura 2000 network), it is difficult to calculate the 
surface areas to evaluate the total contribution of all the 
protected sites. It should be noted that the protected areas 
with the strongest regulatory protection (national parks, 
national and regional nature reserves, and National For- 
est Agency ecological reserves) cover less than 1% of the 
potential niche of this species. 
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Discussion 
Model Performance 


Nine variables, out of the selection of 28 climate vari- 
ables, three topographic variables, and three vegetation 
variables, were found through iteration to contribute 
most effectively to the quality of the models. The tests 
of six models with eight statistical algorithms led to very 
sound results, confirming the performance of the models. 
The consistency between the models and the statistical 
soundness of the modelling are largely due to the good 
spatial coverage of the source data, its geographical pre- 
cision, and the size of the dataset (Araujo and Guisan 
2006). 

Given the exhaustive nature of the occurrence data, we 
feel confident that the predictive model of the ecological 
niche gives a fairly good picture of the potential distribu- 
tion of the species, and can help to map its actual current 
distribution. Other studies have demonstrated that SDMs 
can accurately determine the natural distribution of a spe- 
cies, contributing to the more complete knowledge of its 
current range (Elith and Leathwick 2009). 

In this study, the RF method resulted in the most accu- 
rate SDMs, performing well with all the sets of variables. 
This method uses decision trees based on random group- 
ing of the covariates, modelling both the interactions be- 
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Table 8. Variable importance index generated by the best SDM, sorted by rank of importance. 


Variable Description 

biol4 Precipitation in the driest month 

bio4 Temperature seasonality (std* 100) 

bio9 Mean temperature in the driest quarter 

ndvi Normalized Vegetation Index 

biol5 Precipitation seasonality (coefficient of variation) 
biol0 Mean temperature in the warmest quarter 

biol6 Precipitation in the wettest quarter 

altitude Elevation 

aspect Slope orientation 


tween the variables and their nonlinear relationships, and 
it uses bootstrapping to fit individual trees (Salas et al. 
2017). Cutler (2007) demonstrated that the advantages 
of RF include very high classification accuracy, ability 
to model complex interactions between predictor vari- 
ables, and an algorithm for imputing missing values. 
This gives RF the flexibility to perform several types of 
statistical data analysis, including regression, classifica- 
tion, survival analysis, and unsupervised learning (Cutler 
2007). Rangel and Loyola (2012) also demonstrated that 
machine learning methods such as RF have high statisti- 
cal precision and predictive power for determining the 
species distribution of well-known populations. 

While the models here showed high accuracy, certain 
improvements could be made by integrating additional 
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datasets not included in this study, such as soil maps, the 
distribution of tree species, or detailed vegetation maps 
of France (Leguédois et al. 2011). While land cover 
parameters could contribute to identifying habitat suit- 
ability, their integration in SDMs remains a challenge 
in terms of accuracy and validity (e.g., in terms of how 
current they are). The challenges of incorporating land 
cover data are due to the fact that they are categorical 
variables with limitations imposed by spatial resolution, 
date of production, and thematic classification. Studies 
have shown that continuous remotely sensed predictor 
variables offer many advantages over categorical vari- 
ables and can be used effectively in species distribution 
modelling (Wilson et al. 2013). Models based on biocli- 
matic variables have proven their efficiency in numerous 
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Table 9. The contributions of the different types of protected areas to the conservation of the Ocellated Lizard in France. The surface 
area favorable to the presence of the species (“potential niche” in header of second column) is given in km? according to Model 
M6. The percentage of the protected area relative to the total surface area of the predicted distribution range is shown on the third 


column. 


Surface area under protection within 


Percentage of the potential 


ae the total potential niche (km°) niche* 
Regional Nature Park 8,121 22.68 
Special Area of Conservation (Natura 2000) 5.927 16.55 
Special Protection Areas (Natura 2000) 5.585 15.60 
National Park (park peripheral zone) 792 2.21 
National Nature Reserve 170 0.47 
National Park (park core area) 92 0.26 
Ecological Reserve (National Forest Agency) 43 0.12 
Regional Nature Reserve 20 0.06 


* Total surface area predicted with a probability threshold of 0.7 = 35,805 km”. 


studies, while one study has shown that the addition of 
land cover variables to pure bioclimatic models does not 
necessarily improve the predictive accuracy of the result- 
ing SDM (Thuiller et al. 2004). 


Consistency between Potential and Realized Distribu- 
tions 


The results obtained here showed strong agreement be- 
tween the predictive models and the observed distribu- 
tion of the Ocellated Lizard on a macro-geographic scale, 
giving rise to several conclusions. The first is that climat- 
ic predictors prevail over all the other predictors for this 
species. The current boundaries of the distribution range 
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of the Ocellated Lizard in France are essentially defined 
by climatic factors, which aligns with many studies on 
reptiles (Guisan and Hofer 2003; Santos et al. 2009; Bri- 
to et al. 2011). This suggests an ancient presence of the 
species in France, given the barriers to dispersion such as 
rivers and mountains in its territory, and the time neces- 
sary to colonize the entire potential bioclimatic niche. The 
fragmentation of the populations at the edges of the dis- 
tribution, as well as the historical information regarding 
the loss of populations (Cheylan and Grillet 2005; Doré 
et al. 2015), support this idea and suggest the existence 
in the past of a larger and, above all, a less fragmented 
range. Unfortunately, zooarchaeological information on 
this subject is limited. The species 1s known to have been 
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present in France in the Middle Pleistocene (~700,000 
to 150,000 yr ago), from remains in the Lazaret cave in 
Nice (Bailon 2012), and remains from the Holocene have 
also been found (Mateo 2011). However, the lack of fos- 
sil remains from the last interglacial optimum (between 
125,000 and 11,000 yr ago) does not definitively prove 
the retreat of the species to the Iberian Peninsula dur- 
ing this period. The presence of isolated populations on 
the northern edges of the current distribution, as well as 
its presence in Liguria, would have required overcoming 
major obstacles (the Rhône, Var, and La Roya rivers), 
which points to an ancient occupation of the territory. 
Given these factors, the hypothesis that the species re- 
mained during the interglacial period seems possible, at 
least in the far south of France. 

The strong concurrence between the models and the 
observed distribution also indicates that the process of 
decline in this species is moderate, as all the areas favor- 
able to the species are still occupied, apart from a few 
exceptions. Even at a lower spatial resolution, the bound- 
aries of the distribution range are primarily due to cli- 
matic factors, and only secondarily to ecological factors 
(i.e., presence of favorable habitat). This is particularly 
true at the edges of its distribution range in the valleys 
that open onto the Mediterranean coast (1.e., those of the 
Aude, Rhône, Durance, and Var), where the extent of the 
penetration of the species coincides with the boundary 
of the Mediterranean climate and vegetation, as there 
are no physical obstacles preventing a deeper advance in 
these valleys (Deso et al. 2011, 2015). Notably, the mod- 
el clearly differentiates areas favorable to the Ocellated 
Lizard in zones of rugged terrain. This is particularly the 
case in the region of the Causses (in the southern part of 
Aveyron), which is characterized by limestone plateaus 
that would potentially be favorable to the species but 
where it is not present, and by deep gorges (the valleys of 
the Tarn and the Jonte) where the species has long been 
observed. Rather surprisingly, the model distinguished 
between these two zones (plateaus and gorges) despite 
any notable climatic difference between them. 

On the other hand, several areas not predicted by the 
model have the proven presence of the species (e.g., 
the foothills of the Pyrénées in Ariège, the mountain- 
ous zones of Ardèche, and northern Dordogne). This is 
likely explained by the resolution of the model, which is 
ill-adapted to predicting very small areas within a land- 
scape and climate matrix that is generally unfavorable to 
the species. These known populations live in very small 
micro-habitats (a few dozen ha at most) with unique 
botanical characteristics distinct from the surrounding 
landscapes. An analysis that takes into account a finer 
landscape scale, particularly in terms of vegetation, 
would produce a model with a closer fit. Equally, sub- 
strate characteristics, which were not taken into account 
in the model, play an important role in the presence of 
the Ocellated Lizard when the climatic environment is 
unfavorable. In this case, it seeks out terrain that is rather 
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steep, rocky, or well drained to avoid environments that 
are too wet. 


Which Variables Best Explain the Distribution of the 
Species? 


At the macro-geographic scale, the variables that best ex- 
plain the distribution of the species are related to climate 
and, to a lesser extent, vegetation and topography. This 
indicates the primacy, over all other variables, of a hot 
and dry summer, as well as a strong seasonal contrast; 
two key characteristics of the Mediterranean climate 
(Blondel et al. 2010). The importance of temperature and 
aridity in the summer is certainly due to the reproduction 
requirements of this species. In France, female Ocellated 
Lizards are known to typically lay their eggs at the end 
of May or the beginning of June (Cheylan and Grillet 
2004), and the eggs hatch the third week of September 
or the first week of October (Bischoff et al. 1984; Doré 
et al. 2015). In Provence, this corresponds to an incuba- 
tion period of about 100 days (Cheylan and Grillet 2004). 
Hence, the entire summer period is used for reproduc- 
tion. The late hatching period requires mild temperatures 
at the end of summer and beginning of autumn, allow- 
ing the hatchlings to feed before the hibernation period, 
which begins around 15 November in most of the French 
regions where this species is present (Doré et al. 2015). 

At the local scale, the presence of the species is pri- 
marily influenced by the aridity of the habitat; the Ocel- 
lated Lizard prefers a rocky or well-drained substrate 
that is well exposed to the sun. Dense vegetation cover 
is very unfavorable for this species, as shown in a study 
by Santos and Cheylan (2013) in Provence. In the future, 
gaining a better understanding of the importance of each 
of these habitat variables would be useful, drawing upon 
the resources available on the subject. 


Why is this Species Retreating at the Edges of its 
Distribution Range, in Contrast to Climatic Expecta- 
tions? 


The proven extinction of several Ocellated Lizard popu- 
lations over the last 150 years, mainly on the northern 
border of its distribution range (Cheylan and Grillet 
2005; Grillet et al. 2006) runs counter to what might be 
expected with the warming of the climate, the effects 
of which have been clearly demonstrated on Mediterra- 
nean reptiles in the south of France (Prodon et al. 2017). 
Given its high thermal requirements, this species should 
in fact benefit from climatic warming, particularly at the 
northern edge of its distribution. However, the opposite 
is observed, which suggests the predominance of local 
over global factors. Studies carried out to investigate this 
issue have shown that several local factors explain the 
decline (or even the disappearance) of local populations 
of this species. Those factors include the introduction 
of predators in the case of island populations (Cheylan 
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2016), the disappearance of the European rabbit and the 
resulting changes to the landscape (Grillet et al. 2010), 
the impact of pest control on entomofauna prey (Doré 
et al. 2015), and the abandonment of agricultural land 
and the resulting progression of woodland (Grillet et al. 
2006; Pottier et al. 2017). Thus, the expected effects on 
reptiles of the changes caused by warming in Europe 
(Araújo et al. 2006) are not borne out in the case of the 
Ocellated Lizard. 


Which Zones should be Surveyed in the Future to Im- 
prove Our Knowledge of this Species’ Distribution? 


The recent discovery of a population in Vendée (Cédric 
Baudran, pers. comm. 2018), beyond the known bound- 
ary on the Atlantic coast, shows that new populations 
remain to be discovered, especially at the edges of the 
distribution range. A priority would be to seek confirma- 
tion of the true disappearance of the species in selected 
sites where it was known in the past, based on the pres- 
ence predictions generated by the SDMs. Secondly, an 
attempt to confirm the existence of connections between 
population clusters that are considered to be separated 
would be interesting. This would be particularly useful 
for populations located in the mountainous zones of the 
Alpes-Maritimes (Deso et al. 2015), in the upper Durance 
valley (Deso et al. 2011), and in the Rhône valley (Doré 
et al. 2015), as well as the fragmented populations in the 
Lozere, Aveyron, Tarn, Tarn-et-Garonne, Lot, Dordogne, 
Correze, and Cantal (Geniez and Cheylan 2012; Pottier 
et al. 2017). The coastal populations of the Atlantic cur- 
rently seem rather well-defined (Berroneau 2012); how- 
ever, this does not exclude the possibility of discovering 
new populations there. 


What Conservation Strategy should be Adopted to 
Protect this Species? 


The predictive distribution models generated in this 
study provide interesting leads for defining a conserva- 
tion strategy for this species. First, the current network of 
protected areas in France can be considered to rather sat- 
isfactorily cover the distribution range of the Ocellated 
Lizard and its different population clusters. However, a 
deeper analysis reveals that only a very small proportion 
of the area potentially favorable to the species benefits 
from strong protection regulations. The areas of land 
with the strictest protection (national and regional nature 
reserves, National Forest Agency ecological reserves, 
and national parks) only represent 1.2% of the potential 
niche of this species in France. In terms of national parks, 
the Cévennes National Park clearly bears the most re- 
sponsibility in terms of the conservation of this species, 
followed at some distance by the Calanques National 
Park (respectively, 40 and 95 km? of favorable habitats 
for the species). There are 15 national nature reserves 
with the presence of the species, and in this category of 
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protected area, the reserves of Coussouls de Crau and the 
Maures plain have the largest known populations (re- 
spectively, 74 km? and 52 km? of favorable habitat). Of 
National Forest Agency ecological reserves, 16 include 
land where the lizard is found, with the largest being the 
Maures reserve (18 km? of favorable habitat) and the Pe- 
tit Luberon reserve (16 km? of favorable habitat). Some 
240 Natura 2000 sites have conditions that are potentially 
favorable to the species; 21 protect areas of land that con- 
tribute to the conservation of the species of more than 90 
km?, for a total contribution of about 3,430 km? for this 
category of protected area. 

Given the rather dense network of protected areas, 
both in terms of spatial extent and altitudinal range, a 
strategy based on anticipating climate change (Salas et al. 
2017) is not necessarily the best choice. As stated above, 
this species is in decline at the northern edge of its dis- 
tribution, which runs counter to the expected effects of a 
warming climate (which is predicted for the region in the 
future). Moreover, the refuge habitat for this species (and 
where it originated) is located in the southern half of the 
Iberian Peninsula (Miraldo et al. 2011), so the effects of 
climatic warming are unlikely to be harmful to the spe- 
cies at the northern edge of its distribution. In support 
of this hypothesis, it has been argued that the increase 
in wildfires due to climatic warming will significantly 
increase the density of Ocellated Lizard populations in 
the Mediterranean region, by transforming woodland 
into open landscapes (Santos and Cheylan 2013). A more 
important consideration than climatic warming for the 
conservation of the species is that its spread relies on the 
existence (or not) of favorable environments, and its de- 
mographic capacity to colonize new territories. Unfor- 
tunately, studies of the isolated populations at the edges 
of the distribution range (Grillet et al. 2006; Deso et al. 
2015; Pottier et al. 2017) show that these two parameters 
are rarely present, and that these populations are, in the 
more or less long term, undergoing a process of extinc- 
tion (Salvidio et al. 2004; Cheylan and Grillet 2005). 

From a strategic point of view, therefore, the core of 
the distribution range should be prioritized for conser- 
vation efforts in the long term, without neglecting cer- 
tain peripheral populations in the shorter term (e.g., the 
populations in the valleys of the Durance, Rhône, and 
Var rivers, and the sandy habitats of the Atlantic coast). 
Our SDM-generated maps indicate that the isolated pop- 
ulations of the Atlantic coast, as well as the population 
clusters west of the Massif Central (in the departments 
of Aveyron, Tarn, Tarn-et-Garonne, Lot, and Dordogne), 
offer climatic and topographical conditions that are very 
favorable to this species. These populations, while isolat- 
ed from the Mediterranean population, should be given 
careful attention. They may even harbor specific genetic 
compositions that warrant further consideration. 

The strong dependence of the Ocellated Lizard on 
the European rabbit in soft soils (Grillet et al. 2010) also 
suggests the value of taking concerted conservation mea- 
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sures that equally protect this mammal. As the demo- 
graphic trends of rabbit populations in the Mediterranean 
region are very negative (Ward 2005; Delibes-Mateos 
et al. 2008; Poitevin et al. 2010), this could result in a 
domino effect on Ocellated Lizard populations. 


Conclusions 


This analysis, carried out at the scale of France, reveals 
that the distribution of the Ocellated Lizard is primarily 
conditional on climatic factors, in particular the length of 
the arid summer period. Further study at a smaller scale 
would help to provide a more detailed understanding of 
the ecological preferences of this species. Such a study 
could consider two finer, overlapping spatial scales: the 
Mediterranean coast and the region around Montpellier. 
Focusing on the Mediterranean coastal plains would al- 
low climatic and topographic variables to be separated 
out, at least partially, to better bring to light the roles of 
factors linked to land use. Including the region of Mont- 
pellier would more completely isolate climatic and topo- 
graphic variables, allowing a focus on habitat variables 
directly linked to the ecology of the species: soil type, 
crop or natural vegetation type, level of urbanization, and 
the presence and density of European rabbits. 


Acknowledgements.—We would like to thank all the 
observers and organizations that collected the data used 
in this study: the nature conservation NGOs Cistude-Na- 
ture, the League for the Protection of Birds in Provence- 
Alpes-Céte d’Azur (LPO PACA) and the Drôme (LPO 
Drôme), Nature en Occitanie, Méridionalis, and the pub- 
lic platform for naturalist data, SILENE Provence-Alpes- 
Côte d’ Azur. We would also like to thank Elise Bradbury 
for reviewing the English text. 


Literature Cited 


Allouche O, Tsoar A, Kadmon R. 2006. Assessing the 
accuracy of species distribution models: prevalence, 
kappa, and the true skill statistic (TSS). Journal of Ap- 
plied Ecology 43: 1,223—1,232. 

Anderson RP, Martinez-Meyer H. 2004. Modeling spe- 
cies’ geographic distributions for preliminary conser- 
vation assessments: an implementation with the spiny 
pocket mice (Heteromys) of Ecuador. Biological Con- 
servation 116: 167-179. 

Anderson RP, Lew D, Peterson AT. 2003. Evaluating 
predictive models of species’ distributions: criteria for 
selecting optimal models. Ecological Modelling 162: 
211-232. 

Araujo MB, Segurado P. 2004. An evaluation of methods 
for modelling species distributions. Journal of Bioge- 
ography 31: 1,555-1,568. 

Araujo MB, Guisan A. 2006. Five (or so) challenges for 
species distribution modelling. Journal of Biogeogra- 
phy 33: 1,677-1,688. 


Amphib. Reptile Conserv. 


Araujo MB, Thuiller W, Pearson RG. 2006. Climate 
warming and the decline of amphibians and reptiles 
in Europe. Journal of Biogeography 33(10): 1,712- 
1,728. 

Austin M. 2002. Spatial prediction of species distribu- 
tion: an interface between ecological theory and statis- 
tical modelling. Ecological Modelling 157: 101—118. 

Austin MP, Nicholls AO, Argules CR. 1990. Measure- 
ment of the realized qualitative niche: environmental 
niches of five Eucalyptus species. Ecological Mono- 
graphs 60: 161—177. 

Bailón S. 2012. Données fossiles et mise en place de 
Vherpétofaune actuelle de la France. Pp. 33-39 In: 
Atlas des Amphibiens et Reptiles de France. Coor- 
dinators, Lescure J, Massary de J-C. Biotope, Mèze, 
France and Muséum National d’Histoire Naturelle, 
Paris, France. 272 p. 

Barbet-Massin M, Jiguet F, Albert CH, Thuiller W. 2012. 
Selecting pseudo-absences for species distribution 
models: how, where, and how many? Methods in 
Ecology and Evolution 3(2): 327-338. 

Barbosa AM, Real R, Olivero J, Vargas JM. 2003. Otter 
(Lutra lutra) distribution modeling at two resolution 
scales suited to conservation planning in the Iberian 
Peninsula. Biological Conservation 114(3): 377-387. 

Barry S, Elith J. 2006. Error and uncertainty in habitat 
models. Journal of Applied Ecology 43(3): 413-423. 

Bashfield A, Keim A. 2011. Continent-wide DEM cre- 
ation for the European Union. Pp. 10-15 In: 34th 
International Symposium on Remote Sensing of En- 
vironment. The GEOSS Era: Towards Operational 
Environmental Monitoring. Sydney, Australia. CSI- 
RO, Clayton, Victoria, Australia. 788 p. 

Berroneau M. 2012. Guide Technique de Conservation 
du Lézard Ocellé en Aquitaine. Association Cistude 
Nature, Le Haillan, France. 118 p. 

Bischoff W, Cheylan M, Bohme W. 1984. Lacerta lep- 
ida Daudin 1802 - Perleidechse. Pp. 181—210 In: 
Handbuch der Reptilien und Amphibien Europas: II. 
Editor, Böhme W. Akademische Verlagsgesellschaft, 
Wiesbaden, Germany. 416 p. 

Blondel J, Aronson J, Bodiou JY, Boeuf G. 2010. The 
Mediterranean Region: Biological Diversity in Space 
and Time. Oxford University Press, New York, New 
York, USA. 376 p. 

Bradley A. 1997. The use of the area under the ROC 
curve in the evaluation of machine learning algo- 
rithms. Pattern Recognition 30: 1,145-1,159. 

Breiman L. 2001. Random forests. Machine Learning 
45(1): 5-32. 

Breiman L, Friedman J, Olsen R, Stone C. 1984. Classifi- 
cation and Regression Trees. Wadsworth and Brooks, 
Belmont, California, USA. 358 p. 

Brito JC, Brito-e-Abreu F, Paulo OS, Rosa HD, Crespo 
EG. 1996. Distribution of Schreiber's green lizard 
(Lacerta schreiberi) in Portugal: a predictive model. 
Herpetological Journal 6: 43-47. 


December 2019 | Volume 13 | Number 2 | e213 


Distribution of Timon lepidus in France 


Brito JC, Fahd S, Martinez-Freiria F, Tarroso P, Larbes S, 
Pleguezuelos JM, Santos X. 2011. Climate change and 
peripheral populations: predictions for a relict Medi- 
terranean viper. Acta Herpetologica 6(1): 105-118. 

Brotons L, Thuiller W, Araújo MB, Hirzel AH. 2004. 
Presence-absence versus presence-only modelling 
methods for predicting bird habitat suitability. Ecog- 
raphy 27. 437-448. 

Bucklin DN, Basille M, Benscoter AM, Brandt LA, Maz- 
zotti FJ, Romanach SS, Watling JI. 2015. Comparing 
species distribution models constructed with different 
subsets of environmental predictors. Diversity and 
Distributions 21(1): 23-35. 

Carvalho SB, Brito JC, Crespo EG, Watts ME, Possing- 
ham HP. 2011. Conservation planning under climate 
change: Toward accounting for uncertainty in pre- 
dicted species distributions to increase confidence in 
conservation investments in space and time. Biologi- 
cal Conservation 144(7): 2,020-2,030. 

Cheaib A, Badeau V, Boe J, Chuine I, Delire C, Dufrene 
E, Francois C, Gritti ES, Lesay M, Page C, et al. 
2012. Climate change impacts on tree ranges: Model 
intercomparison facilitates understanding and quan- 
tification of uncertainty. Ecology Letters 15(6): 533— 
544. 

Cheylan M. 2016. The extinction of the Ocellated Lizard 
Timon lepidus lepidus (Daudin, 1802) on the island of 
Porquerolles (Provence, France). Boletin de la Aso- 
ciación Herpetologica Española 27(2): 155—158. 

Cheylan M, Grillet P. 2004. Le Lézard Ocellé. Eveil Na- 
ture. Editions Belin, Paris, France. 95 p. 

Cheylan M, Grillet P. 2005. Statut passé et actuel du Lé- 
zard ocellé (Lacerta lepida, Sauriens, Lacertidae) en 
France. Implication en termes de conservation. Vie et 
Milieu 55(1): 15-30. 

Congedo L, Sallustio L, Munafo M, Ottaviano M, Ton- 
ti D, Marchetti M. 2016. Copernicus high-resolution 
layers for land cover classification in Italy. Journal of 
Maps 12(5): 1,195-1,205. 

Cutler DR, Edwards TC Jr, Beard KH, Cutler A, Hess KT, 
Gibson J, Lawler JJ. 2007. Random forests for classi- 
fication in ecology. Ecology 88(11): 2,783-2,792. 

Delibes-Mateos M, Delibes M, Ferreras P, Villafuerte R. 
2008. Key role of European rabbits in the conserva- 
tion of the western Mediterranean basin hotspot. Con- 
servation Biology 22: 1,106—1,117. 

Deso G, Dusoulier F, Bence S, Cheylan M. 2011. Distri- 
bution du Lézard ocellé Timon lepidus lepidus (Dau- 
din, 1801) dans la Haute vallée de la Durance (Région 
Provence-Alpes-Côte-d’ Azur). Bulletin de la Société 
Herpétologique de France 137: 43—50. 

Deso G, Cevasco J-M, Salvidio S, Ottonella D, Oneto F, 
Cheylan M. 2015. Statut des populations franco-ita- 
liennes de Lézard ocellé Timon lepidus lepidus (Dau- 
din, 1801). Bulletin de la Société Herpétologique de 
France 156: 45-53. 

Doorman CF, Elith J, Bacher S, Buchmann C, Carl G, 


Amphib. Reptile Conserv. 


294 


Carré G, Garcia-Marquez JR, Gruber B, Lafourcade 
B, Leitao PJ, et al. 2012. Collinearity: a review of 
methods to deal with it and a simulation study evalu- 
ating their performance. Ecography 36: 27—46. 

Doré F, Cheylan M, Grillet P. 2015. Le Lézard Ocellé. 
Un Géant sur le Continent Européen. Biotope, Mèze, 
France. 192 p. 

Drake JM, Randin C, Guisan A. 2006. Modelling eco- 
logical niches with support vector machines. Journal 
of Applied Ecology 43: 424—432. 

Elith J, Graham CH, Anderson RP, Dudik M, Ferrier S, 
Guisan A, Hijmans RJ, Huettmann F, Leathwick JR, 
Lehmann A, et al. 2006. Novel methods improve pre- 
diction of species’ distributions from occurrence data. 
Ecography 29(2): 129-151. 

Elith J, Leathwick J. 2007. Predicting species distribu- 
tions from museum and herbarium records using 
multiresponse models fitted with multivariate adap- 
tive regression splines. Diversity and Distribution 13: 
265-275. 

Elith J, Leathwick JR, Hastie T. 2008. A working guide to 
boosted regression trees. Journal of Animal Ecology 
77: 802-813. 

Elith J, Leathwick JR. 2009. Species distribution models: 
ecological explanation and prediction across space 
and time. Annual Review of Ecology, Evolution, and 
Systematics 40: 677-697. 

Engler R, Guisan A, Rechsteiner L. 2004. An improved 
approach for predicting the distribution of rare and 
endangered species from occurrence and pseudo-ab- 
sence data. Journal of Applied Ecology 41(2): 263— 
274. 

Ferreira AF, Quintella BR, Maia C, Mateus CS, Alexan- 
dre CM, Capinha C, Almeida PR. 2013. Influence of 
macrohabitat preferences on the distribution of Euro- 
pean brook and river lampreys: Implications for con- 
servation and management. Biological Conservation 
159: 175-186. 

Ferrier S. 2002. Mapping spatial pattern in biodiversity 
for regional conservation planning: where to from 
here? Systematic Biology 51: 331-363. 

Ferrier S, Watson G, Pearce J, Drielsma M. 2002. Ex- 
tended statistical approaches to modelling spatial pat- 
tern in biodiversity in northeast New South Wales. I. 
Species-level modelling. Biodiversity and Conserva- 
tion 11: 2,275-2,307. 

Fick SE, Hijmans RJ. 2017. WorldClim 2: new 1-km 
spatial resolution climate surfaces for global land 
areas. International Journal of Climatology 37(12): 
4,302-4,315. 

Fielding AH, Bell JF. 1997. A review of methods for the 
assessment of prediction errors in conservation pres- 
ence/absence models. Environmental Conservation 
24(1): 38-49. 

Fourcade Y, Besnard A, Secondi J. 2018. Paintings pre- 
dict the distribution of species, or the challenge of 
selecting environmental predictors and evaluation 


December 2019 | Volume 13 | Number 2 | e213 


Jorcin et al. 


statistics. Global Ecology and Biogeography 27(2): 
245-256. 

Franklin J. 1998. Predicting the distribution of shrub spe- 
cies in southern California from climate and terrain- 
derived variables. Journal of Vegetation Science 9: 
733-748. 

Franklin J. 2009. Mapping Species Distributions, Spa- 
tial Inference, and Prediction. Cambridge University 
Press, Cambridge, United Kingdom. 320 p. 

Geniez P, Cheylan M. 2012. Les Amphibiens et les 
Reptiles du Languedoc-Roussillon et Régions Li- 
mitrophes. Atlas Biogéographique. Biotope, Mèze, 
France, and Muséum National d’Histoire Naturelle, 
Paris, France. 448 p. 

GHRA-LPO Rhone-Alpes (2015). Les Amphibiens et 
Reptiles de Rhône-Alpes. LPO coordination Rhône- 
Alpes, Lyon, France. 448 p. 

Graham CH, Ferrier S, Huettman F, Moritz C, Peterson 
AT. 2004. New developments in museum-based in- 
formatics and applications in biodiversity analysis. 
Trends in Ecology and Evolution 19(9): 497—503. 

Griffith B, Scott JM, Carpenter JW, Reed C. 1989. Trans- 
location as a species conservation tool: status and 
strategy. Science 245: 477-480. 

Grillet P, Cheylan M, Dusoulier F. 2006. Evolution des 
habitats et changements climatique: quelles consé- 
quences pour les populations de lézard ocellé, Lacerta 
lepida (Saurien, Lacertidés) en limite nord de distri- 
bution? Ecologia Mediterranea 32: 63-72. 

Grillet P, Cheylan M, Thirion JM, Doré F, Bonnet X, 
Dauge C, Chollet S, Marchand MA. 2010. Rabbit 
burrows or artificial refuges are a critical habitat com- 
ponent for the threatened lizard Timon lepidus (Sau- 
ria, Lacertidae). Biodiversity and Conservation 19(7): 
2,039-2,051. 

Guisan A, Zimmermann NE. 2000. Predictive habitat 
distribution models in ecology. Ecological Modelling 
135(2-3): 147-186. 

Guisan A, Hofer U. 2003. Predicting reptile distribu- 
tions at the mesoscale: relation to climate and topog- 
raphy. Journal of Biogeography 30(8): 1,233-1,243. 

Guisan A, Thuiller W. 2005. Predicting species distribu- 
tion: offering more than simple habitat models. Ecol- 
ogy Letters 8(9): 993-1 ,009. 

Guisan A, Tingley R, Baumgartner JB, Naujokaitis- 
Lewis I, Sutcliffe PR, Tulloch AIT, Regan TJ, Brotons 
L, McDonald-Madden E, Mantyka-Pringle C, et al. 
2013. Predicting species distributions for conserva- 
tion decisions. Ecology Letters 16(12): 1,424-1,435. 

Harvey PH, Colwell R, Silvertown J, May RM. 2003. 
Null models in ecology. Annual Review of Ecology 
and Systematics 14(1): 189-211. 

Heikkinen RK, Luoto M, Araújo MB. 2006. Methods 
and uncertainties in bioclimatic envelope modelling 
under climate change. Progress in Physical Geogra- 
phy 30(6): 751-777. 

Hirzel A, Guisan A. 2002. Which is the optimal sampling 


Amphib. Reptile Conserv. 


strategy for habitat suitability modelling? Ecological 
Modelling 157: 331-341. 

Huete A, Didan K, van Leeuwen W, Miura T, Glenn E. 
2010. MODIS vegetation indices. Pp. 579—602 In: 
Land Remote Sensing and Global Environmental 
Change. Remote Sensing and Digital Image Process- 
ing, Volume 11. Editors, Ramachandran B, Justice C, 
Abrams M. Springer, New York, New York, USA. 873 


p. 

Humboldt A, Bonpland A. 1805. Essai sur la Géogra- 
phie des Plantes accompagné d'un Tableau Physique 
des Régions Équinoxiales. [Volume] XIII. Levrault, 
Schoell et Compagnie, Paris, France. 155 p. 

Iannella M, Cerasoli F, D’Alessandro P, Console G, 
Biondi M. 2018. Coupling GIS spatial analysis and 
ensemble niche modelling to investigate climate 
change-related threats to the Sicilian pond turtle Emys 
trinacris, an endangered species from the Mediterra- 
nean. PeerJ 6: e4969. 

Jiang G, Sun H, Lang J, Yang L, Li C, Lyet A, Long B, 
Miquelle DG, Zhang C, Aramilev S, et al. 2014. Ef- 
fects of environmental and anthropogenic drivers on 
Amur tiger distribution in northeastern China. Eco- 
logical Research 29(5): 801—813. 

Joly D, Brossard T, Cardot H, Cavailhes J, Hilal M, 
Wavresky P. 2010. Les types de climats en France, 
une construction spatiale. Cybergeo 501: 1-23. 

Karger DN, Conrad O, Böhner J, Kawohl T, Kreft H, 
Soria-Auza RW, Zimmermann NE, Linder HP, Kes- 
sler M. 2017. Climatologies at high resolution for the 
Earth land surface areas. Scientific Data 4: 170122. 

Leguédois S, Party J-P, Dupouey J-L, Gauquelin T, Gé- 
gout J-C, Lecareux C, Badeau V, Probst A. 2011. La 
carte de végétation du CNRS à ľère du numérique. 
European Journal of Geography 559: 1—36. 

Li X, Wang Y. 2012. Applying various algorithms for 
species distribution modelling. /ntegrative Zoology 
8(2): 124-135. 

Liu C, White M, Newell G. 2009. Measuring the ac- 
curacy of species distribution models: a review. Pp. 
4,241-4,247 In: 18th World IMACS/ MODSIM Con- 
gress, Cairns, Australia 13-17 July 2009. Editors, 
Anderssen RS, Braddock RD, Newham LTH. Mod- 
elling and Simulation Society of Australia and New 
Zealand, Canberra, Australia. 

Lyet A, Thuiller W, Cheylan M, Besnard A. 2013. Fine- 
scale regional distribution modelling of rare and 
threatened species: bridging GIS tools and conserva- 
tion in practice. Diversity and Distributions 19(7): 
651-663. 

Manel S, Dias J-M, Ormerdo S. 1999. Comparing dis- 
criminant analysis, neural networks, and logistic re- 
gression for predicting species distributions: a case 
study with a Himalayan river bird. Ecological Model- 
ling 120: 337-347. 

Mateo JA. 2011. Lagarto ocelado - Timon lepidus. In: 
Enciclopedia Virtual de los Vertebrados Españoles. 


December 2019 | Volume 13 | Number 2 | e213 


Distribution of Timon lepidus in France 


Editors, Salvador A, Marco A. Museo Nacional de 
Ciencias Naturales, Madrid, Spain. Available: http:// 
www.vertebradosibericos.org/reptiles/timlep.html 
[Accessed: 14 October 2018]. 

Miraldo A, Faria C, Hewitt GM, Paulo OS, Emerson 
BC. 2013. Genetic analysis of a contact zone between 
two lineages of the Ocellated Lizard (Lacerta lepida 
Daudin 1802) in south-eastern Iberia reveal a steep 
and narrow hybrid zone. Journal of Zoological Sys- 
tematics and Evolutionary Research 51(1): 45—54. 

Moradi S, Ilanloo SS, Kafash A, Yousefi M. 2019. Identi- 
fying high-priority conservation areas for avian biodi- 
versity using species distribution modelling. Ecologi- 
cal Indicators 97: 159-164. 

Muñoz AR, Real R, Barbosa AM, Vargas JM. 2005. 
Modelling the distribution of Bonelli’s eagle in Spain: 
implication for conservation planning. Diversity and 
Distributions 11: 477-486. 

Naimi B, Araújo MB. 2016. Sdm: a reproducible 
and extensible R platform for species distribution 
modelling. Ecography 39: 368-375. 

Naimi B, Hamm NA, Groen TA, Skidmore AK, Toxope- 
us AG. 2014. Where is positional uncertainty a prob- 
lem for species distribution modelling. Ecography 37: 
191-203. 

Nenzen HK, Araujo MB. 2011. Choice of threshold al- 
ters projections of species range shifts under climate 
change. Ecological Modelling 222: 3,346-3,354. 

Poitevin F, Olivier A, Bayle P, Scher O. 2010. Mam- 
miféres de Camargue. Regard du Vivant et Parc Natu- 
rel Régional de Camargue, Castelnau-le-Lez, France. 
232 p. 

Pottier G, Cochard P-O, Barthe L, Cheylan M, Geniez P, 
Defos du Rau P. 2017. État des connaissances sur la 
répartition du Lézard ocellé Timon lepidus (Daudin, 
1802) dans l’ouest de la région Occitanie (ancienne 
région Midi-Pyrénées). Bulletin de la Société Herpé- 
tologique de France 161: 23—56. 

Prodon R, Geniez P, Cheylan M, Devers F, Chuine I, 
Besnard A. 2017. A reversal of the shift towards ear- 
lier spring phenology in several Mediterranean rep- 
tiles and amphibians during the 1998—2013 warming 
slowdown. Global Change Biology 12: 5,481-5,491. 

Pulliam HR. 2000. On the relationship between niche 
and distribution. Ecology Letters 3: 349-361. 

Raes N, ter Steege H. 2007. A null-model for significance 
testing of presence-only species distribution models. 
Ecography 30(5): 121-736. 

Rangel TF, Loyola RD. 2012. Labeling ecological niche 
models. Natureza & Conservação 10(2): 119-126. 
Raxworthy CJ, Martinez-Meyer E, Horning N, Nuss- 
baum RA, Schneider GE, Ortega-Huerta MA, Peter- 
son AT. 2003. Predicting distributions of known and 
unknown reptile species in Madagascar. Nature 426: 

837-841. 

Renwick KM, Curtis C, Kleinhesselink AR, Schlaep- 

fer D, Bradley BA, Aldridge CL, Poulter B, Adler 


Amphib. Reptile Conserv. 


PB. 2018. Multi-model comparison highlights consis- 
tency in predicted effect of warming on a semi-arid 
shrub. Global Change Biology 24: 424-438. 

Ryberg WA, Wolaver BD, Prestridge HL, Labay BJ, 
Pierre JP, Costley RA, Adams S, Bowers BC, Hib- 
bitts TJ. 2017. Habitat modelling and conservation 
of the Western Chicken Turtle (Dierochelys reticu- 
laria miara). Herpetological Conservation and Biol- 
ogy 12(2): 307-320. 

Salas EAL, Seamster VA, Harings NM, Boykin KG, 
Alvarez G, Dixon KW. 2017. Projected future biocli- 
mate-envelope suitability for reptile and amphibian 
species of concern in South Central USA. Herpeto- 
logical Conservation and Biology 12(2): 522—547. 

Salvidio S, Lamagni L, Bombi P, Bologna MA. 2004. 
Distribution, ecology, and conservation status of the 
Ocellated Lizard (Timon lepidus) in Italy (Reptilia, 
Lacertidae). Bollettino di Zoologia 71: 125-134. 

Santos X, Brito JC, Caro J, Abril AJ, Lorenzo M, Sillero 
N, Pleguezuelos JM. 2009. Habitat suitability, threats, 
and conservation of isolated populations of the smooth 
snake (Coronella austriaca) in the southern Iberian 
Peninsula. Biological Conservation 142(2). 344—352. 

Santos X, Cheylan M. 2013. Taxonomic and functional 
response of a Mediterranean reptile community to a 
repeated fire regime. Biological Conservation 168: 
90-98. 

Sillero N. 2011. What does ecological modelling model? 
A proposed classification of ecological niche models 
based on their underlying methods. Ecological Mod- 
elling 222: 1,343—1,346. 

Sillero N, Carretero MA. 2013. Modelling the past and 
future distribution of contracting species. The Iberian 
lizard Podarcis carbonelli (Squamata: Lacertidae) as 
a case study. Zoologischer Anzeiger 252: 289-298. 

Sillero N, Campos J, Bonardi A, Corti C, Creemers R, 
Crochet PA, Crnobrnja Isailović J, Denoël M, Ficetola 
GF, Gonçalves J, et al. 2014. Updated distribution and 
biogeography of amphibians and reptiles of Europe. 
Amphibia-Reptilia 35(1): 1—21. 

Sohrab M, Sayyad SI, Anooshe K, Masoud Y. 2019. Iden- 
tifying high-priority conservation areas for avian bio- 
diversity using species distribution modeling, Eco- 
logical Indicators 97: 159-164. 

Souter NJ, Bull CM, Lethbridge MR, Hutchinson MN. 
2007. Habitat requirements of the endangered Pygmy 
Bluetongue Lizard, Tiliqua adelaidensis. Biological 
Conservation 135(1): 33-45. 

Thuiller W, Araújo MB, Lavorel S. 2004. Do we need 
land-cover data to model species distributions in Eu- 
rope? Journal of Biogeography 31: 353-361. 

Vetaas OR. 2002. Realized and potential climate nich- 
es: a comparison of four Rhododendron tree species. 
Journal of Biogeography 29: 545-554. 

Wan J, Wang C, Yu J, Nie S, Han S, Liu J, Zu Y, Wang 
Q. 2016. Developing conservation strategies for Pinus 
koraiensis and Eleutherococcus senticosus by using 


December 2019 | Volume 13 | Number 2 | e213 


Jorcin et al. 


model-based geographic distributions. Journal of tation structure indices to species distribution models. 
Forestry Research 27(2): 389-400. Biological Conservation 164: 170-176. 

Wang CJ, Wan JZ, Zhang ZX, Zhang GM. 2016. Identi- | Wood SN, Pya N, Safken B. 2016. Smoothing param- 
fying appropriate protected areas for endangered fern eter and model selection for general smooth models. 
species under climate change. SpringerPlus 5(1): 904. Journal of the American Statistical Association 111: 

Ward D. 2005. Reversing Rabbit Decline: One of the 1,548-1,563. 

Biggest Challenges for Nature Conservation in Spain Zomer RJ, Trabucco A, Bossio DA, van Straaten O, Ver- 
and Portugal. Field Reports. IUCN/SSC, Lagomorph chot LV. 2008. Climate change mitigation: A spatial 
Specialist Group, Gland, Switzerland. 54 p. analysis of global land suitability for clean develop- 

Wilson JW, Sexton JO, Jobe RT, Haddad NM. 2013. The ment mechanism afforestation and reforestation. Agri- 

relative contribution of terrain, land cover, and vege- culture Ecosystems and Environment 126: 67—80. 


Laurent Barthe is responsible for biodiversity at the NGO Nature en Occitanie in Toulouse, France. Lau- 
rent has been president of the French Herpetological Society since 2017. Passionate about reptiles since 
childhood, he is particularly interested in snakes and the European Pond Turtle. Two of his job roles are to 
broaden knowledge about reptile and amphibian distribution and promote their conservation. 


Matthieu Berroneau is a herpetologist at the NGO Cistude Nature in Bordeaux, France (http://www. 
cistude.org), where he specializes in the herpetofauna of southwest France, with a particular emphasis on 
conservation, education, and raising awareness. A lifelong interest in amphibians and reptiles also drives 
Matthieu’s work as a professional photographer (http://www.matthieu-berroneau.fr), leading him to travel 
the world to take pictures of a variety of species. 


Marc Cheylan is a lecturer in Conservation Biology at the Ecole Pratique des Hautes Etudes (EPHE), 
an institute within the PSL Research University (Paris Sciences et Lettres). Marc joined the EPHE’s bio- 
geography and ecology research laboratory after a post as the associate curator at the Natural History 
Museum of Aix-en-Provence. He currently works at the French National Centre for Scientific Research 
(CNRS: http://www.cefe.cnrs.fr/fr/), Centre for Functional and Evolutionary Ecology (CEFE), based at the 
University of Montpellier. Author of some 200 publications and four books in the fields of ecology, bio- 
geography, and the conservation of amphibians and reptiles around the Mediterranean, Marc is a member 
of several nature conservation organizations, including the International Union for the Conservation of 
Nature (IUCN), and also sits on various scientific advisory bodies for natural parks and reserves. 


Florian Doré is a naturalist who currently works at the NGO Deux-Sèvres Nature Environnement in 
Deux-Sèvres, western France. Florian leads monitoring surveys and conservation studies on entomofauna 
and herpetofauna, and has contributed to studies of the Ocellated Lizard since 2007. He co-authored the 
National Action Plan for the Ocellated Lizard in France for the NGO OBIOS (Objectifs Biodiversités) 
and co-authored a monograph on the species with Pierre Grillet and Marc Cheylan that was published by 
Biotope Editions. Photo by Marc Cheylan. 


Philippe Geniez is a research engineer in the Vertebrate Biogeography and Ecology lab at the Ecole 
Pratique des Hautes Etudes (EPHE), an institute within the PSL Research University (Paris Sciences 
et Lettres). The lab is part of the French National Centre for Scientific Research (CNRS: http://www. 
cefe.cnrs.fr/fr/), Centre for Functional and Evolutionary Ecology (CEFE). Author of 236 publications, 
Philippe is a specialist in Western Palearctic amphibians and reptiles. His research focuses on biological 
systematics, phylogeny, ecology, and the distribution of plants and animals, particularly amphibians and 
reptiles. 


Pierre Grillet is a naturalist who, since 1995, has been studying the Ocellated Lizard at the edge of its 
distribution, particularly on the island of Oléron, the last island population of the species in France. Pierre 
has written several scientific articles and co-authored two books on the Ocellated Lizard. He regularly 
organizes herpetological training for the French Agency for Biodiversity. Photo by Marc Cheylan. 


Amphib. Reptile Conserv. 297 December 2019 | Volume 13 | Number 2 | e213 


Distribution of Timon lepidus in France 


Pierre Jorcin conducts research and leads projects on biodiversity conservation, with a focus on geospatial 
modelling. Pierre’s areas of interest are species distribution, ecological niche modelling, and wildlife cor- 
ridor mapping. He has been involved in sustainable development programmes in South Asia for 14 years. 
Pierre is currently working on flora and fauna database management for ecological ranking and environ- 
mental impact assessment studies in southern France. 


Benjamin Kabouche is the director of the environmental NGO LPO PACA (Ligue pour la Protection 
des Oiseaux, a member of BirdLife International) in the region of Provence-Alpes-Céte d’ Azur in France 
(https://paca.lpo.fr/protection). One of Benjamin’s roles is to coordinate studies and nature conservation 
programs. He has contributed to several naturalist publications on the subjects of terrestrial wildlife and 
biogeography. 


Alexandre Movia works for the environmental NGO LPO Drôme (Ligue pour la Protection des Oiseaux), 
where he is an ecological corridor specialist. Alexandre acts as an advisor on herpetology to the French 
department of the Drôme, and recently conducted a study on the distribution of the Ocellated Lizard in 
this region. 


Babak Naimi is a researcher at the University of Helsinki, Finland, with a research focus on modelling 
species distribution and biodiversity under climate change and land use change scenarios. Babak is inter- 
ested in developing a quantitative understanding of ecosystem dynamics (e.g., through remote sensing and 
geoinformatics tools) and uncovering the complexities behind ecosystem behavior. 


Gilles Pottier has been a professional field herpetologist for some 20 years, and currently works mainly 
in the Pyrenees and the Massif Central for the NGO Nature en Occitanie. A member of the French Herpe- 
tological Society, Gilles runs training events in herpetology in southwest France and has written numer- 
ous papers and books about the local herpetofauna, including a work on the reptiles of the Pyrenees (Les 
Reptiles des Pyrénées), published in 2016 by the French Natural History Museum. Photo by J.P. Vacher. 


Jean-Marc Thirion is an ecologist and the director of the NGO OBIOS (Objectifs Biodiversités) in south- 
west France. Jean-Marc leads conservation projects to protect natural areas and conducts population moni- 
toring of amphibians and reptiles to promote conservation initiatives. He has participated in many natural- 
ist surveys of flora and fauna. 


Amphib. Reptile Conserv. 298 December 2019 | Volume 13 | Number 2 | e213 


